function plotSparsityLP(n_x,n_u,n_w,n_p_opt,n_tp,n_eq_Ceq, ...
               n_eq_Cineq,n_eq_Heq,n_eq_Hineq,n_e,n_cp,n_eq_F0, ...
               n_eq_F, dF0_n_cols,dF_dp_n_cols,dF_ddx_dx_du_dw_n_cols, ...
               dCeq_dp_n_cols,dCeq_ddx_dx_du_dw_n_cols,dCeq_ddx_p_dx_p_du_p_dw_p_n_cols, ... 			  
           dCineq_dp_n_cols,dCineq_ddx_dx_du_dw_n_cols,dCineq_ddx_p_dx_p_du_p_dw_p_n_cols, ...			  
           dHeq_dp_n_cols,dHeq_ddx_p_dx_p_du_p_dw_p_n_cols, ...			  
           dHineq_dp_n_cols,dHineq_ddx_p_dx_p_du_p_dw_p_n_cols, ...              
           ind_dg,ind_dh,fig_nbr)

%    Copyright (C) 2009 Modelon AB
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, version 3 of the License.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.


n_rows = n_eq_F0 + ...
    	 n_eq_F*n_e*n_cp + ...
         n_x*n_e + ...
         n_x*n_e*n_cp + ...
	 n_u + ...
         (2*n_x  + n_u + n_w)*n_tp + ...
         n_eq_Ceq*(n_e*n_cp + 1) + ...
         n_eq_Heq

n_cols = n_p_opt + ...
         (2*n_x + n_u + n_w)*(n_e*n_cp + 1) + ...
         n_x*n_e + ...
         (2*n_x + n_u + n_w)*n_tp       

col_F0 = [0.8 0.8 0.8];
col_F = [0.8 0.8 0.8];
col_el_junc = [0.3 0.3 0.3];
col_cp = [0.5 0.5 0.5];


figure(fig_nbr);
clf
p_F0 = patch([0 dF0_n_cols+1 dF0_n_cols+1 0], ...
      [0 0 n_eq_F0+1 n_eq_F0+1],col_F0);  
hold on
set(p_F0,'EdgeAlpha',0);

plot([0 0],[0 n_rows+1],'Color',col_el_junc)
plot([0 n_cols + 1],[0 0],'Color',col_el_junc)

for i=1:n_e*n_cp,
  p_F = patch([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*(i-1) ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i ...
	dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i ...
	       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*(i-1)], ...
       [n_eq_F0 + 1 + n_eq_F*(i-1) ...
        n_eq_F0 + 1 + n_eq_F*(i-1) ...
        n_eq_F0 + 1 + n_eq_F*(i) ...
        n_eq_F0 + 1 + n_eq_F*(i)],col_F);
  set(p_F,'EdgeAlpha',0);
end;




% plot([dF0_n_cols + 1 dF0_n_cols + 1], [0 n_rows+1])

for i=0:n_e,
        
  if (i~=n_e), 
    for j=1:n_cp-1,
      plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp + ...
	  dF_ddx_dx_du_dw_n_cols*j ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp + ...
	  dF_ddx_dx_du_dw_n_cols*j], ...
       [0 n_rows+1], 'Color', col_cp)
    end
  end
  
  plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp], ...
       [0 n_rows+1], 'Color', col_el_junc)
  
		       
end;


%for i=0:n_e*n_cp,
%  plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*i ...
%		      n_eq_F0 + 1 + n_eq_F*i], ...
%       'Color',col_cp);
%end


for i=1:n_e,
  plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
	n_x*i ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
       n_x*i], ...
       [0 n_rows+1], 'Color', col_el_junc)  
end


  plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
	n_x*n_e + (n_x+n_x+n_w+n_u)*n_tp ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
       n_x*n_e + (n_x+n_x+n_w+n_u)*n_tp ], ...
       [0 n_rows+1], 'Color', col_el_junc)  

  plot([0,n_cols+1], [n_eq_F0 + 1 ...
		      n_eq_F0 + 1], ...
      'Color',col_el_junc);


  plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp ...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp], ...
      'Color',col_el_junc);

  
    plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
			n_x*n_e  ...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
		       n_x*n_e ], ...
      'Color',col_el_junc);

       plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
			n_x*n_e+n_u...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
		       n_x*n_e+n_u], ...
      'Color',col_el_junc);

    
    plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
			n_x*n_e + n_x*n_e*n_cp + n_u...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
		       n_x*n_e + n_x*n_e*n_cp + n_u], ...
      'Color',col_el_junc);

    plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
			n_x*n_e + n_x*n_e*n_cp + n_u + (n_x+n_x+n_w+n_u)*n_tp...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
		       n_x*n_e + n_x*n_e*n_cp + n_u+ (n_x+n_x+n_w+n_u)*n_tp], ...
      'Color',col_el_junc);

      plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
			n_x*n_e + n_x*n_e*n_cp + n_u + (n_x+n_x+n_w+n_u)*n_tp + n_eq_Ceq*(n_e*n_cp + 1)...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
		       n_x*n_e + n_x*n_e*n_cp + n_u + (n_x+n_x+n_w+n_u)*n_tp+ n_eq_Ceq*(n_e*n_cp + 1)], ...
      'Color',col_el_junc);

        plot([0,n_cols+1], [n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
			n_x*n_e + n_x*n_e*n_cp + n_u + (n_x+n_x+n_w+n_u)*n_tp + n_eq_Ceq*(n_e*n_cp + 1) + ...
            n_eq_Heq ...
		      n_eq_F0 + 1 + n_eq_F*n_e*n_cp + ...
		       n_x*n_e + n_x*n_e*n_cp + n_u + (n_x+n_x+n_w+n_u)*n_tp+ n_eq_Ceq*(n_e*n_cp + 1) + ...
               n_eq_Heq], ...
      'Color',col_el_junc);

  
plot(ind_dh(:,3)+0.5,ind_dh(:,2)+0.5,'x');
set(gca,'YDir','reverse');                                     

title('Sparsity matrix of the equality constraint Jacobian')

figure(fig_nbr+1)
clf

plot(ind_dg(:,3)+0.5,ind_dg(:,2)+0.5,'x');
set(gca,'YDir','reverse');                                     

hold on

title('Sparsity matrix of the inequality constraint Jacobian')

n_rows = n_eq_Cineq*(n_e*n_cp+1) + n_eq_Hineq;


plot([0 0],[0 n_rows+1],'Color',col_el_junc)
plot([0 n_cols + 1],[0 0],'Color',col_el_junc)
plot([0 n_cols + 1],[n_eq_Cineq*(n_e*n_cp+1)+1 n_eq_Cineq*(n_e*n_cp+1)+1],'Color',col_el_junc)
plot([0 n_cols + 1],[n_rows+1 n_rows+1],'Color',col_el_junc)

for i=0:n_e,
        
  if (i~=n_e), 
    for j=1:n_cp-1,
      plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp + ...
	  dF_ddx_dx_du_dw_n_cols*j ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp + ...
	  dF_ddx_dx_du_dw_n_cols*j], ...
       [0 n_rows+1], 'Color', col_cp)
    end
  end
  
  plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*i*n_cp], ...
       [0 n_rows+1], 'Color', col_el_junc)
  
		       
end;

for i=1:n_e,
  plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
	n_x*i ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
       n_x*i], ...
       [0 n_rows+1], 'Color', col_el_junc)  
end


  plot([dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
	n_x*n_e + (n_x+n_x+n_w+n_u)*n_tp ...
       dF0_n_cols + 1 + dF_ddx_dx_du_dw_n_cols*n_e*n_cp + ...
       n_x*n_e + (n_x+n_x+n_w+n_u)*n_tp ], ...
       [0 n_rows+1], 'Color', col_el_junc)  


